library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(knitr)
library(jagsUI)
## Warning: package 'jagsUI' was built under R version 4.4.3
library(DiagrammeR)
## Warning: package 'DiagrammeR' was built under R version 4.4.3
set.seed(1019)
Seed is set to 1019 for reproducibility.
Manual data, pen-wise, inclusive of multiple broods.
dflongforcount <- data.frame(
dose=factor(c(rep("cont",6),rep("low",7),rep("med",6),rep("high",8)),
levels=c("cont","low","med","high")),
dosenum=c(rep(0.000121,6),rep(0.121,7),rep(0.45,6),rep(2.04,8)),
pen=c(31,32,33,34,35,35,36,37,38,38,39,40,40,41,42,43,43,44,45,46,46,46,47,48,48,48,50),
clutch=c(1,1,1,1,1,2,1,1,1,2,1,1,2,1,1,1,2,1,1,1,2,3,1,1,2,3,1),
numhatched=c(3,2,2,1,0,0,4,1,0,3,3,0,0,3,0,0,0,1,1,0,0,0,2,0,0,0,2)
)
dflongforcount$failure <- ifelse(dflongforcount$numhatched==0,1,0)
jags.data <- list(y=dflongforcount$numhatched,
conc=dflongforcount$dosenum,
n=nrow(dflongforcount))
cat(file="model_AMKE_poisson_eggs hatched by normal dose.txt","
model{
#priors and linear models
alpha ~ dnorm(log(3),100) # mean and precision in jags precision is inverse of the variance, so 1/sigma^2 and in this case suggests the priors are mean=0, sd=0.1 (exp(0.1) on the real scale) (sqrt(1/100))
beta ~ dnorm(0,0.001)
#likelihood for the poisson GLM
for(i in 1:n){
y[i] ~ dpois(lambda[i]) #stochastic
log(lambda[i]) <- alpha + beta * conc[i]
}
mean.exp.count <- exp(alpha)
}
")
inits <- function(){list(alpha=rnorm(1,0,1),beta=rnorm(1,0,1))}
parameters <- c("alpha","beta","mean.exp.count","lambda")
ni <- 50000
nb <- 10000
nc <- 3
nt <- 10
na <- 1000
out1 <- jags(jags.data, inits, parameters, "model_AMKE_poisson_eggs hatched by normal dose.txt", n.iter=ni, n.burnin=nb, n.chains=nc, n.thin=nt, n.adapt=na, parallel=T)
##
## Processing function input.......
##
## Done.
##
## Beginning parallel processing using 3 cores. Console output will be suppressed.
##
## Parallel processing completed.
##
## Calculating statistics.......
##
## Done.
test <- data.frame(out1$sims.list[1:2])
#write_csv(test, "out1posterioralphaandbeta.csv")
\[ \alpha \sim N(3,0.1)\\ \beta \sim N(0,31.62)\\ y \sim Poisson(\lambda)\\ log(\lambda) = alpha + \beta \sum PFAS \]
contm <- out1$mean$lambda[1] #2.575
lowm <- out1$mean$lambda[7] #2.264
medm <- out1$mean$lambda[14] #1.605
highm <- out1$mean$lambda[20] #0.343
contl <- out1$q2.5$lambda[1] #2.164
lowl <- out1$q2.5$lambda[7] #1.906
medl <- out1$q2.5$lambda[14] #1.181
highl <- out1$q2.5$lambda[20] #0.081
contu <- out1$q97.5$lambda[1] #3.041
lowu <- out1$q97.5$lambda[7] #2.66
medu <- out1$q97.5$lambda[14] #2.034
highu <- out1$q97.5$lambda[20] #0.812
### probability of having a certain number of viable eggs
counteggdf <- data.frame(
dose=c(0.000121, 0.121, 0.45, 2.04),
mean=c(contm, lowm, medm, highm),
lower=c(contl, lowl, medl, highl),
upper=c(contu, lowu, medu, highu)
)
effectzerodf <- data.frame(
dose=c(0.000121, 0.121, 0.45, 2.04),
mean=c(dpois(0, lambda=contm),dpois(0, lambda=lowm),dpois(0, lambda=medm),dpois(0, lambda=highm)),
lower=c(dpois(0, lambda=contl),dpois(0, lambda=lowl),dpois(0, lambda=medl),dpois(0, lambda=highl)),
upper=c(dpois(0, lambda=contu),dpois(0, lambda=lowu),dpois(0, lambda=medu),dpois(0, lambda=highu)),
hatchedeggs=rep(0,4)
)
effect3df <- data.frame(
dose=c(0.000121, 0.121, 0.45, 2.04),
mean=c(dpois(3, lambda=contm),dpois(3, lambda=lowm),dpois(3, lambda=medm),dpois(3, lambda=highm)),
lower=c(dpois(3, lambda=contl),dpois(3, lambda=lowl),dpois(3, lambda=medl),dpois(3, lambda=highl)),
upper=c(dpois(3, lambda=contu),dpois(3, lambda=lowu),dpois(3, lambda=medu),dpois(3, lambda=highu)),
hatchedeggs=rep(3,4)
)
effect2df <- data.frame(
dose=c(0.000121, 0.121, 0.45, 2.04),
mean=c(dpois(2, lambda=contm),dpois(2, lambda=lowm),dpois(2, lambda=medm),dpois(2, lambda=highm)),
lower=c(dpois(2, lambda=contl),dpois(2, lambda=lowl),dpois(2, lambda=medl),dpois(2, lambda=highl)),
upper=c(dpois(2, lambda=contu),dpois(2, lambda=lowu),dpois(2, lambda=medu),dpois(2, lambda=highu)),
hatchedeggs=rep(2,4)
)
effect5df <- data.frame(
dose=c(0.000121, 0.121, 0.45, 2.04),
mean=c(dpois(5, lambda=contm),dpois(5, lambda=lowm),dpois(5, lambda=medm),dpois(5, lambda=highm)),
lower=c(dpois(5, lambda=contl),dpois(5, lambda=lowl),dpois(5, lambda=medl),dpois(5, lambda=highl)),
upper=c(dpois(5, lambda=contu),dpois(5, lambda=lowu),dpois(5, lambda=medu),dpois(5, lambda=highu)),
hatchedeggs=rep(5,4)
)
effect1df <- data.frame(
dose=c(0.000121, 0.121, 0.45, 2.04),
mean=c(dpois(1, lambda=contm),dpois(1, lambda=lowm),dpois(1, lambda=medm),dpois(1, lambda=highm)),
lower=c(dpois(1, lambda=contl),dpois(1, lambda=lowl),dpois(1, lambda=medl),dpois(1, lambda=highl)),
upper=c(dpois(1, lambda=contu),dpois(1, lambda=lowu),dpois(1, lambda=medu),dpois(1, lambda=highu)),
hatchedeggs=rep(1,4)
)
effect4df <- data.frame(
dose=c(0.000121, 0.121, 0.45, 2.04),
mean=c(dpois(4, lambda=contm),dpois(4, lambda=lowm),dpois(4, lambda=medm),dpois(4, lambda=highm)),
lower=c(dpois(4, lambda=contl),dpois(4, lambda=lowl),dpois(4, lambda=medl),dpois(4, lambda=highl)),
upper=c(dpois(4, lambda=contu),dpois(4, lambda=lowu),dpois(4, lambda=medu),dpois(4, lambda=highu)),
hatchedeggs=rep(4,4)
)
bigeggprobdf <- rbind(
effectzerodf,
effect1df,
effect2df,
effect3df,
effect4df,
effect5df
)
Probability of having a certain number of hatched eggs at each dose
with confidence intervals.
Note that upper and lower are strange because some dose response are
decreasing and some are increasing depending on the number of hatched
eggs.
bigeggprobdf |>
#filter(hatchedeggs==0)|>
kable(, digits=5)
| dose | mean | lower | upper | hatchedeggs |
|---|---|---|---|---|
| 0.00012 | 0.07615 | 0.11519 | 0.04752 | 0 |
| 0.12100 | 0.10421 | 0.14819 | 0.06999 | 0 |
| 0.45000 | 0.20227 | 0.30665 | 0.13041 | 0 |
| 2.04000 | 0.71318 | 0.92599 | 0.44442 | 0 |
| 0.00012 | 0.19609 | 0.24894 | 0.14477 | 1 |
| 0.12100 | 0.23566 | 0.28293 | 0.18614 | 1 |
| 0.45000 | 0.32326 | 0.36247 | 0.26566 | 1 |
| 2.04000 | 0.24107 | 0.07120 | 0.36042 | 1 |
| 0.00012 | 0.25247 | 0.26901 | 0.22053 | 2 |
| 0.12100 | 0.26645 | 0.27010 | 0.24750 | 2 |
| 0.45000 | 0.25831 | 0.21423 | 0.27058 | 2 |
| 2.04000 | 0.04074 | 0.00274 | 0.14615 | 2 |
| 0.00012 | 0.21671 | 0.19379 | 0.22396 | 3 |
| 0.12100 | 0.20084 | 0.17190 | 0.21940 | 3 |
| 0.45000 | 0.13761 | 0.08441 | 0.18373 | 3 |
| 2.04000 | 0.00459 | 0.00007 | 0.03951 | 3 |
| 0.00012 | 0.13951 | 0.10470 | 0.17058 | 4 |
| 0.12100 | 0.11354 | 0.08205 | 0.14586 | 4 |
| 0.45000 | 0.05498 | 0.02495 | 0.09357 | 4 |
| 2.04000 | 0.00039 | 0.00000 | 0.00801 | 4 |
| 0.00012 | 0.07185 | 0.04526 | 0.10394 | 5 |
| 0.12100 | 0.05135 | 0.03133 | 0.07758 | 5 |
| 0.45000 | 0.01757 | 0.00590 | 0.03812 | 5 |
| 2.04000 | 0.00003 | 0.00000 | 0.00130 | 5 |
Using the posterior draws, find the doses bracketing 95% probability of observing that amount of nest failure.
dpois(0,lambda=contm)
## [1] 0.07614831
dpois(0,lambda=contm)+0.2
## [1] 0.2761483
i.e. “safe” concentration
### do prediction interpolation for having zero
alt.conc.pred <- seq(from=0.0001,to=3, length.out=1000)
pred2 <- array(NA,dim=c(length(alt.conc.pred),12000))
for(i in 1:12000){
pred2[,i] <- exp(out1$sims.list$alpha[i]+out1$sims.list$beta[i]*alt.conc.pred)
}
pm2 <- apply(pred2,1, mean)
pm2l <- apply(pred2,1, quantile, probs=0.975)
pm2u <- apply(pred2,1, quantile, probs=0.025)
prob <- apply(pred2, 1, function(x) mean((x-1)<0))
dprob0 <- apply(pred2,1,function(x) mean(dpois(0,lambda=x)))
dprob0min <- apply(pred2,1,function(x) quantile(dpois(0,lambda=x), probs=0.05))
dprob0max <- apply(pred2,1,function(x) quantile(dpois(0,lambda=x), probs=0.95))
#plot(dprob0~alt.conc.pred, ylim=c(0,1))
#points(dprob0min~alt.conc.pred, col="darkgray")
#points(dprob0max~alt.conc.pred, col="darkgray")
mintwentythreshold <- max(alt.conc.pred[dprob0min<0.28])
maxtwentythreshold <- max(alt.conc.pred[dprob0max<0.28])
meantwentythreshold <- max(alt.conc.pred[dprob0<0.28])
mintwentythreshold
## [1] 1.090154
meantwentythreshold
## [1] 0.648727
maxtwentythreshold
## [1] 0.432518
par(mai=c(2,1.2,0.1,0.3))
plot(jitter(numhatched, 0.5)~jitter(as.numeric(factor(dose))), data=dflongforcount, pch=3, ylab="Number of hatched eggs", xlab=expression(Adult~daily~dose~mg/kg~Sigma*PFAS),
bty="n", las=1, axes=F, cex.lab=1.5)
axis(1, at=c(1,2,3,4), labels=c("0.000121", "0.121", "0.45", "2.04"), font=1, cex.axis=1.5)
axis(2, at=c(0,1,2,3,4), labels=c(0,1,2,3,4), font=1, las=1, cex.axis=1.5)
arrows(as.numeric(factor(counteggdf$dose)),counteggdf$lower,as.numeric(factor(counteggdf$dose)),counteggdf$upper, lwd=2, length=0.1, angle=90, code=3, col=c("darkgreen","blue","darkorange","red"))
points(mean~as.numeric(factor(dose)), data=counteggdf, type="l", col="darkgray", lwd=2)
points(mean~as.numeric(factor(dose)), data=counteggdf, type="p", col=c("darkgreen","blue","darkorange","red"), pch=16, cex=1.5)
points(c(1.1,1.1),c(2.95,3.05),col="purple", pch=2)
axis(1, at=c(1,2,3,4), labels=c("0.003", "3.35", "8.34", "65.84"), font=1, cex.axis=1.5, line=5)
mtext(expression(Mean~egg~mg/kg~Sigma*PFAS), side=1, line=8, font=1, cex=1.5)
par(mai=c(1,1,0.1,0.1))
plot(pm2~alt.conc.pred, log="", ylim=c(0,3), type="l", bty="n", xlim=c(-0.25,3),
ylab="Mean number of hatched eggs", xlab=expression(paste("Adult daily dose ", Sigma,"PFAS mg/kg")), las=1,
font=1, font.lab=1, pch=16, cex=2, lwd=2, cex.lab=1.5, cex.axis=1.5)
points(pm2l~alt.conc.pred, type="l", lty=2, col="darkgray", lwd=2)
points(pm2u~alt.conc.pred, type="l", lty=2, col="darkgray", lwd=2)
text(-0.35,2.575,"2.575", pos=4, col="darkgreen")
text(-0.35,1.264456,"1.264", pos=4, col="red")
segments(0,1.264456,1.5,1.264456, col="red")
segments(0.6757532,1.264456,0.6757532,0.5, col="black")
text(0.6757532,0.5,"0.68 mg/kg", pos=1, col="black", xpd=NA)
segments(0.4054919,1.264456,0.4054919,1.064456, col="black")
#text(0.432518,0.48,"0.43 mg/kg\nP(0 eggs)<0.28", pos=3, col="black")
text(0.4054919,1.064456,"0.41 mg/kg", pos=1, col="black")
segments(1.237296,1.264456,1.237296,1.464456, col="black")
#text(1.090154,0.18,"1.09 mg/kg\nP(0 eggs)<0.28", pos=1, col="black", xpd=NA)
text(1.237296,1.464456,"1.24 mg/kg", pos=3, col="black", xpd=NA)
max(alt.conc.pred[pm2>1.264456])
## [1] 0.6757532
max(alt.conc.pred[pm2l>1.264456])
## [1] 1.237296
max(alt.conc.pred[pm2u>1.264456])
## [1] 0.4054919
par(mai=c(1.2,1.2,0.1,0.75))
plot(dprob0~alt.conc.pred, log="", ylim=c(0,1), type="l", bty="n", xlim=c(-0.25,3),
ylab="Probability of having 0 hatched eggs", xlab=expression(paste("Adult daily dose ", Sigma,"PFAS mg/kg")), las=1,
font=1, font.lab=1, pch=16, cex=2, lwd=2, cex.lab=1.5, cex.axis=1.5)
points(dprob0min~alt.conc.pred, type="l", lty=2, col="darkgray", lwd=2)
points(dprob0max~alt.conc.pred, type="l", lty=2, col="darkgray", lwd=2)
abline(h=c(0,1))
text(-0.35,dpois(0,lambda=contm),"0.08", pos=4, col="darkgreen")
text(-0.35,dpois(0,lambda=contm)+0.2,"0.28", pos=4, col="red")
segments(0,dpois(0,lambda=contm)+0.2,1.5,dpois(0,lambda=contm)+0.2, col="red")
segments(0.432518,0.28,0.432518,0.48, col="black")
#text(0.432518,0.48,"0.43 mg/kg\nP(0 eggs)<0.28", pos=3, col="black")
text(0.432518,0.48,"0.43 mg/kg", pos=3, col="black")
segments(1.090154,0.28,1.090154,0.18, col="black")
#text(1.090154,0.18,"1.09 mg/kg\nP(0 eggs)<0.28", pos=1, col="black", xpd=NA)
text(1.090154,0.18,"1.09 mg/kg", pos=1, col="black", xpd=NA)
segments(0.648727,0.28,0.648727,0.08, col="black")
text(0.648727,0.08,"0.65 mg/kg", pos=1, col="black", xpd=NA)
jags.data2 <- list(y=dflongforcount$failure,
conc=dflongforcount$dosenum,
n=nrow(dflongforcount))
cat(file="model_AMKE_bernoulli_eggs hatched by normal dose.txt","
model{
#priors and linear models
alpha <- logit(mean.theta)
mean.theta ~ dunif(0,0.5)
beta ~ dnorm(0,0.001)
#likelihood for the bernoulli GLM
for(i in 1: n){
y[i] ~ dbern(theta[i]) #stochastic
logit(theta[i]) <- alpha + beta * conc[i]
}
}
")
inits2 <- function(){list(mean.theta=runif(1,0, 0.5),beta=rnorm(1,0,1))}
parameters2 <- c("alpha","beta","mean.theta","theta")
ni <- 50000
nb <- 10000
nc <- 3
nt <- 10
na <- 1000
out2 <- jags(jags.data2, inits2, parameters2, "model_AMKE_bernoulli_eggs hatched by normal dose.txt", n.iter=ni, n.burnin=nb, n.chains=nc, n.thin=nt, n.adapt=na, parallel=T)
##
## Processing function input.......
##
## Done.
##
## Beginning parallel processing using 3 cores. Console output will be suppressed.
##
## Parallel processing completed.
##
## Calculating statistics.......
##
## Done.
#jagsUI::traceplot(out2)
test2 <- data.frame(out2$sims.list[1:2])
#write_csv(test, "out1posterioralphaandbeta.csv")
\[ \alpha \sim \text{Uniform}(0,0.5)\\ \beta \sim N(0,31.62)\\ y \sim \text{Bernoulli}(\theta)\\ logit(\theta) = alpha + \beta \sum PFAS \]
contb <- out2$mean$theta[1]
lowb <- out2$mean$theta[7]
medb <- out2$mean$theta[14]
highb <- out2$mean$theta[21]
contbl <- out2$q2.5$theta[1]
lowbl <- out2$q2.5$theta[7]
medbl <- out2$q2.5$theta[14]
highbl <- out2$q2.5$theta[21]
contbu <- out2$q97.5$theta[1]
lowbu <- out2$q97.5$theta[7]
medbu <- out2$q97.5$theta[14]
highbu <- out2$q97.5$theta[21]
### probability of a failed nest (having zero hatched eggs)
failuredf <- data.frame(
dose=c(0.000121, 0.121, 0.45, 2.04),
mean=c(contb, lowb, medb, highb),
lower=c(contbl, lowbl, medbl, highbl),
upper=c(contbu, lowbu, medbu, highbu)
)
# plot(jitter(dflongforcount$failure, 0.2)~jitter(as.numeric(factor(dflongforcount$dosenum)),0.5), axes=F, pch=3, ylab="Per nest probability of failure", xlab="Diet concentration (mg/kg-d sumPFAS)", )
# lines(as.numeric(factor(dflongforcount$dosenum)), out2$mean$theta, col="blue", lwd=3)
# lines(as.numeric(factor(dflongforcount$dosenum)), out2$q2.5$theta, lwd=3, lty=2)
# lines(as.numeric(factor(dflongforcount$dosenum)), out2$q97.5$theta, lwd=3, lty=2)
# axis(1,at=unique(as.numeric(factor(dflongforcount$dosenum))), labels=sort(unique(round(dflongforcount$dosenum,2))))
# axis(2,las=1)
par(mai=c(2,1.2,0.1,0.3))
plot(jitter(failure, 0.5)~jitter(as.numeric(factor(dose))), data=dflongforcount, pch=3, ylab="Probability of nest failure", xlab=expression(Adult~daily~dose~mg/kg~Sigma*PFAS),
bty="n", las=1, axes=F, cex.lab=1.5)
axis(1, at=c(1,2,3,4), labels=c("0.000121", "0.121", "0.45", "2.04"), font=1, cex.axis=1.5)
axis(2, at=c(0,0.5,1), labels=c(0,0.5,1), font=1, las=1, cex.axis=1.5)
arrows(as.numeric(factor(failuredf$dose)),failuredf$lower,as.numeric(factor(failuredf$dose)),failuredf$upper, lwd=2, length=0.1, angle=90, code=3, col=c("darkgreen","blue","darkorange","red"))
points(mean~as.numeric(factor(dose)), data=failuredf, type="l", col="darkgray", lwd=2)
points(mean~as.numeric(factor(dose)), data=failuredf, type="p", col=c("darkgreen","blue","darkorange","red"), pch=16, cex=1.5)
points(c(0.95,1.05),c(0,0),col="purple", pch=2)
axis(1, at=c(1,2,3,4), labels=c("0.003", "3.35", "8.34", "65.84"), font=1, cex.axis=1.5, line=5)
mtext(expression(Mean~egg~mg/kg~Sigma*PFAS), side=1, line=8, font=1, cex=1.5)
Using the posterior draws, find the doses bracketing 95% probability of observing that amount of nest failure.
dbinom(1,1,prob=contb)
## [1] 0.3577885
dbinom(1,1,prob=contb)+0.2
## [1] 0.5577885
i.e. “safe” concentration
### do prediction interpolation for having zero
alt.conc.pred <- seq(from=0.0001,to=3, length.out=1000)
pred3 <- array(NA,dim=c(length(alt.conc.pred),12000))
for(i in 1:12000){
pred3[,i] <- plogis(out2$sims.list$alpha[i]+out2$sims.list$beta[i]*alt.conc.pred)
}
pm3 <- apply(pred3,1, mean)
pm3l <- apply(pred3,1, quantile, probs=0.975)
pm3u <- apply(pred3,1, quantile, probs=0.025)
prob2 <- apply(pred3, 1, function(x) mean((x-1)<0))
dprob02 <- apply(pred3,1,function(x) mean(dbinom(1,1,prob=x)))
dprob02min <- apply(pred3,1,function(x) quantile(dbinom(1,1,prob=x), probs=0.05))
dprob02max <- apply(pred3,1,function(x) quantile(dbinom(1,1,prob=x), probs=0.95))
#plot(dprob0~alt.conc.pred, ylim=c(0,1))
#points(dprob0min~alt.conc.pred, col="darkgray")
#points(dprob0max~alt.conc.pred, col="darkgray")
mintwentythreshold2 <- max(alt.conc.pred[dprob02min<0.5579766])
maxtwentythreshold2 <- max(alt.conc.pred[dprob02max<0.5579766])
meantwentythreshold2 <- max(alt.conc.pred[dprob02<0.5579766])
mintwentythreshold2
## [1] 2.915919
meantwentythreshold2
## [1] 0.8919622
maxtwentythreshold2
## [1] 0.3574454
par(mai=c(1.2,1.2,0.1,0.75))
plot(dprob02~alt.conc.pred, log="", ylim=c(0,1), type="l", bty="n", xlim=c(-0.25,3),
ylab="Probability of nest failure", xlab=expression(paste("Adult daily dose ", Sigma,"PFAS mg/kg")), las=1,
font=1, font.lab=1, pch=16, cex=2, lwd=2, cex.lab=1.5, cex.axis=1.5)
points(dprob02min~alt.conc.pred, type="l", lty=2, col="darkgray", lwd=2)
points(dprob02max~alt.conc.pred, type="l", lty=2, col="darkgray", lwd=2)
abline(h=c(0,1))
text(-0.35,dbinom(1,1,prob=contb),"0.36", pos=4, col="darkgreen")
text(-0.35,dbinom(1,1,prob=contb)+0.2,"0.56", pos=4, col="red")
segments(0,dbinom(1,1,prob=contb)+0.2,3,dbinom(1,1,prob=contb)+0.2, col="red")
segments(0.3514396,0.56,0.3514396,0.66, col="black")
text(0.3514396,0.66,"0.35 mg/kg", pos=3, col="black")
segments(3,0.56,3,0.26, col="black")
text(2.915919,0.26,"2.92 mg/kg", pos=1, col="black", xpd=NA)
segments(0.8919622,0.56,0.8919622,0.36, col="black")
text(0.8919622,0.36,"0.89 mg/kg", pos=1, col="black", xpd=NA)
mermaid("
graph RL
A(Adult);
J(Juvenile);
I(Immigration);
A--Fecundity-->J;
A--Adult Survival-->A;
J--Juvenile Survival-->A;
I-->A;
",
width=500, height=500)
Survival and immigration rates from McClure et al. 2021. McClure CJW et al. 2021. Demography of a widespread raptor across disparate regions. Ibis. 163(2):658–670 https://doi.org/10.1111/ibi.12916
juvsurvdist <- c(
rnorm(1000, mean=0.13, sd=0.06),
rnorm(1000, mean=0.15, sd=0.05),
rnorm(1000, mean=0.09, sd=0.03),
rnorm(1000, mean=0.08, sd=0.06)
)
juvsurvdist <- juvsurvdist[!juvsurvdist<0]
par(mfrow=c(1,2),mai=c(0.8,0.8,0.1,0.1))
hist(juvsurvdist, main="")
plot(density(juvsurvdist), main="",bty="n")
par(mfrow=c(1,1),mai=c(1,1,0.5,0.5))
summary(juvsurvdist)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0003585 0.0756261 0.1106364 0.1159136 0.1525845 0.3493614
sd(juvsurvdist)
## [1] 0.05639346
Across the populations.
sd(c(0.13, 0.15, 0.09, 0.08))
## [1] 0.03304038
adtsurvdist <- c(
rnorm(1000, mean=0.50, sd=0.03),
rnorm(1000, mean=0.48, sd=0.03),
rnorm(1000, mean=0.55, sd=0.03),
rnorm(1000, mean=0.56, sd=0.08)
)
#adtsurvdist <- adtsurvdist[adtsurvdist>0]
par(mfrow=c(1,2),mai=c(0.8,0.8,0.1,0.1))
hist(adtsurvdist, main="")
plot(density(adtsurvdist), main="",bty="n")
par(mfrow=c(1,1),mai=c(1,1,0.5,0.5))
summary(adtsurvdist)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.3272 0.4826 0.5151 0.5229 0.5548 0.7859
sd(adtsurvdist)
## [1] 0.05953504
sd(c(0.50, 0.48, 0.55, 0.56))
## [1] 0.0386221
immidist <- c(
rnorm(1000, mean=0.21, sd=0.06),
rnorm(1000, mean=0.21, sd=0.05),
rnorm(1000, mean=0.15, sd=0.04),
rnorm(1000, mean=0.20, sd=0.09)
)
immidist <- immidist[!immidist<0]
par(mfrow=c(1,2),mai=c(0.8,0.8,0.1,0.1))
hist(immidist, main="")
plot(density(immidist), main="",bty="n")
par(mfrow=c(1,1),mai=c(1,1,0.5,0.5))
summary(immidist)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.001566 0.148561 0.191600 0.194103 0.236840 0.474021
sd(immidist)
## [1] 0.06687488
sd(c(0.21, 0.21, 0.15, 0.20))
## [1] 0.02872281
fecunddist <- c(
rnorm(1000, mean=1.55, sd=0.13),
rnorm(1000, mean=1.36, sd=0.05),
rnorm(1000, mean=1.38, sd=0.06),
rnorm(1000, mean=1.32, sd=0.13)
)
#fecunddist <- fecunddist[!fecunddist<0]
par(mfrow=c(1,2),mai=c(0.8,0.8,0.1,0.1))
hist(fecunddist, main="")
plot(density(fecunddist), main="",bty="n")
par(mfrow=c(1,1),mai=c(1,1,0.5,0.5))
summary(fecunddist)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.9304 1.3272 1.3820 1.4032 1.4625 1.9546
sd(fecunddist)
## [1] 0.1328875
sd(c(1.55, 1.36, 1.38, 1.32))
## [1] 0.1014479
See https://github.com/mikemeredith/IPM_code/blob/main/IPM_03/IPM_03.3.5.R
Code has been modified so “temporal” stochasticity is replaced with “population” stochasticity based on the 4 populations in Mcclure et al. 2021.
# https://github.com/kenkellner/IPMbook/blob/main/R/stopifnot.R
stopifnotProbability <- function(arg, allowNA=FALSE) {
name <- deparse(substitute(arg))
if(allowNA && all(is.na(arg))) { # An all-NA vector is logical, but ok.
# do nothing
} else {
if(!allowNA && any(is.na(arg)))
stop("Argument '", name, "' must not contain NA or NaN.", call.=FALSE)
if(!is.numeric(arg))
stop("Argument '", name, "' must be numeric.", call.=FALSE)
if(any(arg < 0 | arg > 1, na.rm=TRUE)) {
if(allowNA) {
stop("Argument '", name, "' must be a probability between 0 and 1, or NA.", call.=FALSE)
} else {
stop("Argument '", name, "' must be a probability between 0 and 1.", call.=FALSE)
}
}
}
}
stopifNegative <- function(arg, allowNA=FALSE, allowZero=TRUE) {
name <- deparse(substitute(arg))
if(allowNA && all(is.na(arg))) { # An all-NA vector is logical, but ok.
# do nothing
} else {
if(!allowNA && any(is.na(arg)))
stop("Argument '", name, "' must not contain NA or NaN.", call.=FALSE)
if(!is.numeric(arg))
stop("Argument '", name, "' must be numeric.", call.=FALSE)
if(allowZero) {
if(any(arg < 0, na.rm=TRUE)) {
if(allowNA) {
stop("Argument '", name, "' must be non-negative, or NA.", call.=FALSE)
} else {
stop("Argument '", name, "' must be non-negative.", call.=FALSE)
}
}
} else {
if(any(arg <= 0, na.rm=TRUE)) {
if(allowNA) {
stop("Argument '", name, "' must be greater than 0, or NA.", call.=FALSE)
} else {
stop("Argument '", name, "' must be greater than 0.", call.=FALSE)
}
}
}
}
}
# https://github.com/kenkellner/IPMbook/blob/main/R/stopifnot.R
# https://github.com/kenkellner/IPMbook/blob/main/R/BetaDist.R
getBeta2Par <- function(mean, sd) {
stopifnotProbability(mean, allowNA=FALSE)
stopifNegative(sd, allowNA=FALSE, allowZero=FALSE)
nu <- mean * (1-mean) / sd^2 - 1
if(any(nu <= 0)) {
warning("sd is too large; some shape parameters will be NA.", call.=FALSE)
nu[nu <= 0] <- NA
}
alpha <- mean * nu
beta <- (1-mean) * nu
cbind(shape1=alpha, shape2=beta)
}
dbeta2 <- function(x, mean, sd) {
shapes <- getBeta2Par(mean,sd)
return(dbeta(x, shapes[,1], shapes[,2]))
}
pbeta2 <- function(q, mean, sd, lower.tail=TRUE, log.p=FALSE) {
shapes <- getBeta2Par(mean,sd)
return(pbeta(q, shapes[,1], shapes[,2], lower.tail=lower.tail, log.p=log.p))
}
qbeta2 <- function(p, mean, sd, lower.tail=TRUE, log.p=FALSE) {
shapes <- getBeta2Par(mean,sd)
return(qbeta(p, shapes[,1], shapes[,2], lower.tail=lower.tail, log.p=log.p))
}
rbeta2 <- function(n, mean, sd) {
shapes <- getBeta2Par(mean,sd)
return(rbeta(n, shapes[,1], shapes[,2]))
}
# https://github.com/kenkellner/IPMbook/blob/main/R/BetaDist.R
# https://github.com/mikemeredith/IPM_code/blob/main/IPM_03/IPM_03.3.5.R
# 3.3 Classical analysis of a matrix population model
# ===================================================
# 3.3.5 Analysis of a matrix population model with different sources of
# stochasticity and parameter uncertainty
# ---------------------------------------------------------------------
# Define mean, measurement error and temporal variability of the demographic parameters
mean.sj <- 0.1163615 # Mean value of juv. survival
sd.sj.e <- 0.0554588 # Uncertainty of mean juv. survival as SD on natural scale
sd.sj.t <- plogis(0.03304038) # population variability of juv. survival as SD on logit scale # Andrew, this might have to be plogis(0.03304038), but not sure, check the text.
mean.sa <- 0.5222 # Mean value of ad. survival
sd.sa.e <- 0.05769693 # Uncertainty of mean ad. survival as SD on natural scale
sd.sa.t <- plogis(0.0386221) # Temporal variability of ad. survival as SD on logit scale
mean.f1 <- 1.4035 # Mean value of productivity of 1y females
sd.f1.e <- 0.1326032 # Uncertainty of mean productivity as SD on natural scale
sd.f1.t <- 0.1014479 # Temporal variability of productivity as SD on natural scale
mean.fa <- 1.4035 # Mean value of productivity of adult females
sd.fa.e <- 0.1326032 # Uncertainty of mean productivity as SD on natural scale
sd.fa.t <- 0.101 # Temporal variability of productivity as SD on natural scale
mean.imm <- 0.192587 # Mean value of adult immigration
sd.imm.e <- 0.06674809 # Uncertainty of mean adult immigration as SD on natural scale
sd.imm.t <- 0.02872281 # population variability of adult immigration as SD on natrual scale
# Define the number of years with predictions and the Monte Carlo setting
T <- 10 # Number of years (projection time frame)
nsim <- 10000 # Number of replicate populations simulated
# Define population matrix and initial stage-specific population sizes
N <- array(NA, dim=c(2, T+1, nsim))
N[,1,] <- c(3,1)
r <- matrix(NA, nrow=T, ncol=nsim)
alive <- matrix(NA, nrow=T, ncol=nsim)
mean.r <- numeric(nsim)
# Project population
for (s in 1:nsim){ # Loop over replicate populations
#if(s %% 100 == 0) {cat(paste("*** Simrep", s, "***\n")) } # Counter
# Generate a mean of the demographic rates (subject to measurement error)
msj <- rbeta2(1, mean.sj, sd.sj.e)
msa <- rbeta2(1, mean.sa, sd.sa.e)
mf1 <- rnorm(1, mean.f1, sd.f1.e)
mfa <- rnorm(1, mean.fa, sd.fa.e)
mia <- rnorm(1, mean.imm, sd.imm.e)
# Generate annual demographic rates (subject to temporal variability)
sj <- plogis(rnorm(T, qlogis(msj), sd.sj.t))
sa <- plogis(rnorm(T, qlogis(msa), sd.sa.t))
f1 <- pmax(0, rnorm(T, mf1, sd.f1.t)) # Avoids negative values
fa <- pmax(0, rnorm(T, mfa, sd.fa.t)) # Ditto
ia <- pmax(0, rnorm(T, mia, sd.imm.t)) # Ditto
# Project population (include demographic stochasticity)
for (t in 1:T){ # Loop over years
N[1,t+1,s] <- rpois(1, sj[t] * (f1[t] * N[1,t,s] + fa[t] * N[2,t,s]))
N[2,t+1,s] <- rbinom(1, sum(N[,t,s]), sa[t])+rpois(1,ia)
if (sum(N[,t+1,s]) == 0) break
#r[t,s] <- log(sum(N[,t+1,s])) - log(sum(N[,t,s]))
alive[t,s] <- t
} #t
#mean.r[s] <- mean(r[min(alive[,s], na.rm=TRUE):max(alive[,s], na.rm=TRUE),s])
}
#mean(mean.r)
#sd(mean.r)
not.extinct <- which(!is.na(alive[T,]))
mean(mean.r[not.extinct])
## [1] 0
sd(mean.r[not.extinct])
## [1] 0
# Extinction probability (after T years)
sum(is.na(alive[T,])) / nsim
## [1] 0.8031
rowSums(is.na(alive)) / nsim
## [1] 0.0342 0.1379 0.2686 0.3941 0.5065 0.5914 0.6652 0.7232 0.7672 0.8031
df <- data.frame(
totalcount=c(colSums(N[,,1:nsim])),
time=rep(c(0:T),nsim),
simnumber=gl(nsim, T+1)
)
# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# boxplot(c(colSums(N[,,1:nsim]))~factor(rep(c(0:T),nsim)))
bxpobj<-boxplot(c(colSums(N[,,1:nsim]))~factor(rep(c(0:T),nsim)), plot=F)
# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# points(bxpobj$stats[3,]~c(0:T), type="l", col="red", lwd=2)
#
# par(mai=c(1,1,0.1,1))
# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# points(bxpobj$stats[3,]~c(0:T), type="b", col="blue", lwd=2, pch=16)
# par(new=T)
# plot(c(rowSums(is.na(alive))/nsim)~c(1:T), type="b", col="red", bty="n", xlab="", ylab="", axes=F, ylim=c(0,1), xlim=c(0,T), lwd=2, pch=16)
# axis(side=4, at=pretty(range(c(0,1))))
# mtext("probability of extinction", side=4, line=3)
par(mai=c(1,1,0.1,1))
plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3, xlab="years of simulation", ylab="estimated population size")
points(bxpobj$stats[3,]~c(0:T), type="b", col="blue", lwd=2, pch=16)
par(new=T)
plot(c(rowSums(is.na(alive))/nsim)~c(1:T), type="b", col="red", bty="n", xlab="", ylab="", axes=F, ylim=c(0,1), xlim=c(0,T), lwd=2, pch=16)
axis(side=4, at=pretty(range(c(0,1))))
mtext("proportion of simulations reached N=0", side=4, line=3, col="red")
ggplot()+
geom_line(data=df, aes(x=time, y=totalcount, group=simnumber), linewidth=1, color="gray", alpha=0.8) +
geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
#geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
geom_label(aes(x=c(1:T), y=0, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
theme_classic(base_size=22) +
labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")
## Warning: Removed 40883 rows containing missing values or values outside the scale range
## (`geom_line()`).
ggplot()+
geom_jitter(data=df, aes(x=time, y=totalcount, group=simnumber), shape=".", width=0.3, height=0.3) +
geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
#geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
geom_label(aes(x=c(1:T), y=-0.5, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
theme_classic(base_size=22)+
labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")
## Warning: Removed 40883 rows containing missing values or values outside the scale range
## (`geom_point()`).
# ggplot()+
# geom_hex(data=df, aes(x=time, y=totalcount, group=simnumber), bins=10) +
# geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
# #geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
# geom_label(aes(x=c(1:T), y=-0.5, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
# scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
# theme_classic(base_size=22)+
# labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")
# ~~~~ Plot graph with population growth rates (Fig. 3.15) ~~~~
# op <- par(mar=c(4, 4, 2, 1), las=1, cex=1.1, "mfrow")
# layout(matrix(c(1, 1, 2, 3), 2, 2, byrow=TRUE), widths=c(1, 1), heights=c(1, 1), TRUE)
# plot(r[,1], type="l", lwd=0.5, ylab="Annual population growth rate", xlab="Time",
# ylim=range(r[which(!is.na(alive))]), col="lightgrey", axes=FALSE)
# axis(1); axis(2)
# for (s in 2:nsim){
# lines(r[!is.na(alive[,s]),s], lwd=0.5, col="lightgrey")
# }
# lines(apply(r, 1, mean, na.rm=TRUE), lwd=1.5)
# mtext("A", at=0, cex=1.5)
# a <- hist(mean.r, nclass=50, col="dodgerblue", main="",
# xlab="Population growth rate", axes=FALSE) #xlim=c(-2.5, 0.1),
# axis(1)
# axis(2, at=c(0, 5000, 10000, 15000, 20000, 25000, 30000), labels=c(0, 5, 10, 15, 20, 25, 30))
# mtext("B", at=a$mids[1], cex=1.5)
# a <- hist(mean.r[not.extinct], nclass=25, col="dodgerblue", main="",
# xlab="Population growth rate", axes=FALSE)
# axis(1)
# axis(2, at=c(0, 2000, 4000, 6000, 8000, 10000), labels=c(0, 2, 4, 6, 8, 10))
# mtext("C", at=a$mids[1], cex=1.5)
# par(op)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# https://github.com/mikemeredith/IPM_code/blob/main/IPM_03/IPM_03.3.5.R
# https://github.com/mikemeredith/IPM_code/blob/main/IPM_03/IPM_03.3.5.R
# 3.3 Classical analysis of a matrix population model
# ===================================================
# 3.3.5 Analysis of a matrix population model with different sources of
# stochasticity and parameter uncertainty
# ---------------------------------------------------------------------
# Define mean, measurement error and temporal variability of the demographic parameters
mean.sj <- 0.1163615 # Mean value of juv. survival
sd.sj.e <- 0.0554588 # Uncertainty of mean juv. survival as SD on natural scale
sd.sj.t <- plogis(0.03304038) # population variability of juv. survival as SD on logit scale # Andrew, this might have to be plogis(0.03304038), but not sure, check the text.
mean.sa <- 0.5222 # Mean value of ad. survival
sd.sa.e <- 0.05769693 # Uncertainty of mean ad. survival as SD on natural scale
sd.sa.t <- plogis(0.0386221) # Temporal variability of ad. survival as SD on logit scale
mean.f1 <- 2.575 # Mean value of productivity of 1y females
sd.f1.e <- 0.1326032 # Uncertainty of mean productivity as SD on natural scale
sd.f1.t <- 0.1014479 # Temporal variability of productivity as SD on natural scale
mean.fa <- 2.575 # Mean value of productivity of adult females
sd.fa.e <- 0.1326032 # Uncertainty of mean productivity as SD on natural scale
sd.fa.t <- 0.101 # Temporal variability of productivity as SD on natural scale
mean.imm <- 0.192587 # Mean value of adult immigration
sd.imm.e <- 0.06674809 # Uncertainty of mean adult immigration as SD on natural scale
sd.imm.t <- 0.02872281 # population variability of adult immigration as SD on natrual scale
# Define the number of years with predictions and the Monte Carlo setting
T <- 10 # Number of years (projection time frame)
nsim <- 10000 # Number of replicate populations simulated
# Define population matrix and initial stage-specific population sizes
N <- array(NA, dim=c(2, T+1, nsim))
N[,1,] <- c(3,1)
r <- matrix(NA, nrow=T, ncol=nsim)
alive <- matrix(NA, nrow=T, ncol=nsim)
mean.r <- numeric(nsim)
# Project population
for (s in 1:nsim){ # Loop over replicate populations
#if(s %% 100 == 0) {cat(paste("*** Simrep", s, "***\n")) } # Counter
# Generate a mean of the demographic rates (subject to measurement error)
msj <- rbeta2(1, mean.sj, sd.sj.e)
msa <- rbeta2(1, mean.sa, sd.sa.e)
mf1 <- rnorm(1, mean.f1, sd.f1.e)
mfa <- rnorm(1, mean.fa, sd.fa.e)
mia <- rnorm(1, mean.imm, sd.imm.e)
# Generate annual demographic rates (subject to temporal variability)
sj <- plogis(rnorm(T, qlogis(msj), sd.sj.t))
sa <- plogis(rnorm(T, qlogis(msa), sd.sa.t))
f1 <- pmax(0, rnorm(T, mf1, sd.f1.t)) # Avoids negative values
fa <- pmax(0, rnorm(T, mfa, sd.fa.t)) # Ditto
ia <- pmax(0, rnorm(T, mia, sd.imm.t)) # Ditto
# Project population (include demographic stochasticity)
for (t in 1:T){ # Loop over years
N[1,t+1,s] <- rpois(1, sj[t] * (f1[t] * N[1,t,s] + fa[t] * N[2,t,s]))
N[2,t+1,s] <- rbinom(1, sum(N[,t,s]), sa[t])+rpois(1,ia)
if (sum(N[,t+1,s]) == 0) break
#r[t,s] <- log(sum(N[,t+1,s])) - log(sum(N[,t,s]))
alive[t,s] <- t
} #t
#mean.r[s] <- mean(r[min(alive[,s], na.rm=TRUE):max(alive[,s], na.rm=TRUE),s])
}
#mean(mean.r)
#sd(mean.r)
not.extinct <- which(!is.na(alive[T,]))
mean(mean.r[not.extinct])
## [1] 0
sd(mean.r[not.extinct])
## [1] 0
# Extinction probability (after T years)
sum(is.na(alive[T,])) / nsim
## [1] 0.5943
rowSums(is.na(alive)) / nsim
## [1] 0.0254 0.0902 0.1752 0.2637 0.3443 0.4097 0.4697 0.5164 0.5589 0.5943
df <- data.frame(
totalcount=c(colSums(N[,,1:nsim])),
time=rep(c(0:T),nsim),
simnumber=gl(nsim, T+1)
)
# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# boxplot(c(colSums(N[,,1:nsim]))~factor(rep(c(0:T),nsim)))
bxpobj<-boxplot(c(colSums(N[,,1:nsim]))~factor(rep(c(0:T),nsim)), plot=F)
# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# points(bxpobj$stats[3,]~c(0:T), type="l", col="red", lwd=2)
#
# par(mai=c(1,1,0.1,1))
# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# points(bxpobj$stats[3,]~c(0:T), type="b", col="blue", lwd=2, pch=16)
# par(new=T)
# plot(c(rowSums(is.na(alive))/nsim)~c(1:T), type="b", col="red", bty="n", xlab="", ylab="", axes=F, ylim=c(0,1), xlim=c(0,T), lwd=2, pch=16)
# axis(side=4, at=pretty(range(c(0,1))))
# mtext("probability of extinction", side=4, line=3)
par(mai=c(1,1,0.1,1))
plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3, xlab="years of simulation", ylab="estimated population size")
points(bxpobj$stats[3,]~c(0:T), type="b", col="blue", lwd=2, pch=16)
par(new=T)
plot(c(rowSums(is.na(alive))/nsim)~c(1:T), type="b", col="red", bty="n", xlab="", ylab="", axes=F, ylim=c(0,1), xlim=c(0,T), lwd=2, pch=16)
axis(side=4, at=pretty(range(c(0,1))))
mtext("proportion of simulations reached N=0", side=4, line=3, col="red")
ggplot()+
geom_line(data=df, aes(x=time, y=totalcount, group=simnumber), linewidth=1, color="gray", alpha=0.8) +
geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
#geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
geom_label(aes(x=c(1:T), y=0, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
theme_classic(base_size=22) +
labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")
## Warning: Removed 28535 rows containing missing values or values outside the scale range
## (`geom_line()`).
ggplot()+
geom_jitter(data=df, aes(x=time, y=totalcount, group=simnumber), shape=".", width=0.3, height=0.3) +
geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
#geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
geom_label(aes(x=c(1:T), y=-0.5, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
theme_classic(base_size=22)+
labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")
## Warning: Removed 28535 rows containing missing values or values outside the scale range
## (`geom_point()`).
# ggplot()+
# geom_hex(data=df, aes(x=time, y=totalcount, group=simnumber), bins=10) +
# geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
# #geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
# geom_label(aes(x=c(1:T), y=-0.5, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
# scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
# theme_classic(base_size=22)+
# labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")
# ~~~~ Plot graph with population growth rates (Fig. 3.15) ~~~~
# op <- par(mar=c(4, 4, 2, 1), las=1, cex=1.1, "mfrow")
# layout(matrix(c(1, 1, 2, 3), 2, 2, byrow=TRUE), widths=c(1, 1), heights=c(1, 1), TRUE)
# plot(r[,1], type="l", lwd=0.5, ylab="Annual population growth rate", xlab="Time",
# ylim=range(r[which(!is.na(alive))]), col="lightgrey", axes=FALSE)
# axis(1); axis(2)
# for (s in 2:nsim){
# lines(r[!is.na(alive[,s]),s], lwd=0.5, col="lightgrey")
# }
# lines(apply(r, 1, mean, na.rm=TRUE), lwd=1.5)
# mtext("A", at=0, cex=1.5)
# a <- hist(mean.r, nclass=50, col="dodgerblue", main="",
# xlab="Population growth rate", axes=FALSE) #xlim=c(-2.5, 0.1),
# axis(1)
# axis(2, at=c(0, 5000, 10000, 15000, 20000, 25000, 30000), labels=c(0, 5, 10, 15, 20, 25, 30))
# mtext("B", at=a$mids[1], cex=1.5)
# a <- hist(mean.r[not.extinct], nclass=25, col="dodgerblue", main="",
# xlab="Population growth rate", axes=FALSE)
# axis(1)
# axis(2, at=c(0, 2000, 4000, 6000, 8000, 10000), labels=c(0, 2, 4, 6, 8, 10))
# mtext("C", at=a$mids[1], cex=1.5)
# par(op)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# https://github.com/mikemeredith/IPM_code/blob/main/IPM_03/IPM_03.3.5.R
# https://github.com/mikemeredith/IPM_code/blob/main/IPM_03/IPM_03.3.5.R
# 3.3 Classical analysis of a matrix population model
# ===================================================
# 3.3.5 Analysis of a matrix population model with different sources of
# stochasticity and parameter uncertainty
# ---------------------------------------------------------------------
# Define mean, measurement error and temporal variability of the demographic parameters
mean.sj <- 0.1163615 # Mean value of juv. survival
sd.sj.e <- 0.0554588 # Uncertainty of mean juv. survival as SD on natural scale
sd.sj.t <- plogis(0.03304038) # population variability of juv. survival as SD on logit scale # Andrew, this might have to be plogis(0.03304038), but not sure, check the text.
mean.sa <- 0.5222 # Mean value of ad. survival
sd.sa.e <- 0.05769693 # Uncertainty of mean ad. survival as SD on natural scale
sd.sa.t <- plogis(0.0386221) # Temporal variability of ad. survival as SD on logit scale
mean.f1 <- 2.264 # Mean value of productivity of 1y females
sd.f1.e <- 0.1326032 # Uncertainty of mean productivity as SD on natural scale
sd.f1.t <- 0.1014479 # Temporal variability of productivity as SD on natural scale
mean.fa <- 2.264 # Mean value of productivity of adult females
sd.fa.e <- 0.1326032 # Uncertainty of mean productivity as SD on natural scale
sd.fa.t <- 0.101 # Temporal variability of productivity as SD on natural scale
mean.imm <- 0.192587 # Mean value of adult immigration
sd.imm.e <- 0.06674809 # Uncertainty of mean adult immigration as SD on natural scale
sd.imm.t <- 0.02872281 # population variability of adult immigration as SD on natrual scale
# Define the number of years with predictions and the Monte Carlo setting
T <- 10 # Number of years (projection time frame)
nsim <- 10000 # Number of replicate populations simulated
# Define population matrix and initial stage-specific population sizes
N <- array(NA, dim=c(2, T+1, nsim))
N[,1,] <- c(3,1)
r <- matrix(NA, nrow=T, ncol=nsim)
alive <- matrix(NA, nrow=T, ncol=nsim)
mean.r <- numeric(nsim)
# Project population
for (s in 1:nsim){ # Loop over replicate populations
#if(s %% 100 == 0) {cat(paste("*** Simrep", s, "***\n")) } # Counter
# Generate a mean of the demographic rates (subject to measurement error)
msj <- rbeta2(1, mean.sj, sd.sj.e)
msa <- rbeta2(1, mean.sa, sd.sa.e)
mf1 <- rnorm(1, mean.f1, sd.f1.e)
mfa <- rnorm(1, mean.fa, sd.fa.e)
mia <- rnorm(1, mean.imm, sd.imm.e)
# Generate annual demographic rates (subject to temporal variability)
sj <- plogis(rnorm(T, qlogis(msj), sd.sj.t))
sa <- plogis(rnorm(T, qlogis(msa), sd.sa.t))
f1 <- pmax(0, rnorm(T, mf1, sd.f1.t)) # Avoids negative values
fa <- pmax(0, rnorm(T, mfa, sd.fa.t)) # Ditto
ia <- pmax(0, rnorm(T, mia, sd.imm.t)) # Ditto
# Project population (include demographic stochasticity)
for (t in 1:T){ # Loop over years
N[1,t+1,s] <- rpois(1, sj[t] * (f1[t] * N[1,t,s] + fa[t] * N[2,t,s]))
N[2,t+1,s] <- rbinom(1, sum(N[,t,s]), sa[t])+rpois(1,ia)
if (sum(N[,t+1,s]) == 0) break
#r[t,s] <- log(sum(N[,t+1,s])) - log(sum(N[,t,s]))
alive[t,s] <- t
} #t
#mean.r[s] <- mean(r[min(alive[,s], na.rm=TRUE):max(alive[,s], na.rm=TRUE),s])
}
#mean(mean.r)
#sd(mean.r)
not.extinct <- which(!is.na(alive[T,]))
mean(mean.r[not.extinct])
## [1] 0
sd(mean.r[not.extinct])
## [1] 0
# Extinction probability (after T years)
sum(is.na(alive[T,])) / nsim
## [1] 0.6536
rowSums(is.na(alive)) / nsim
## [1] 0.0270 0.1089 0.2069 0.2994 0.3825 0.4560 0.5194 0.5701 0.6143 0.6536
df <- data.frame(
totalcount=c(colSums(N[,,1:nsim])),
time=rep(c(0:T),nsim),
simnumber=gl(nsim, T+1)
)
# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# boxplot(c(colSums(N[,,1:nsim]))~factor(rep(c(0:T),nsim)))
bxpobj<-boxplot(c(colSums(N[,,1:nsim]))~factor(rep(c(0:T),nsim)), plot=F)
# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# points(bxpobj$stats[3,]~c(0:T), type="l", col="red", lwd=2)
#
# par(mai=c(1,1,0.1,1))
# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# points(bxpobj$stats[3,]~c(0:T), type="b", col="blue", lwd=2, pch=16)
# par(new=T)
# plot(c(rowSums(is.na(alive))/nsim)~c(1:T), type="b", col="red", bty="n", xlab="", ylab="", axes=F, ylim=c(0,1), xlim=c(0,T), lwd=2, pch=16)
# axis(side=4, at=pretty(range(c(0,1))))
# mtext("probability of extinction", side=4, line=3)
par(mai=c(1,1,0.1,1))
plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3, xlab="years of simulation", ylab="estimated population size")
points(bxpobj$stats[3,]~c(0:T), type="b", col="blue", lwd=2, pch=16)
par(new=T)
plot(c(rowSums(is.na(alive))/nsim)~c(1:T), type="b", col="red", bty="n", xlab="", ylab="", axes=F, ylim=c(0,1), xlim=c(0,T), lwd=2, pch=16)
axis(side=4, at=pretty(range(c(0,1))))
mtext("proportion of simulations reached N=0", side=4, line=3, col="red")
ggplot()+
geom_line(data=df, aes(x=time, y=totalcount, group=simnumber), linewidth=1, color="gray", alpha=0.8) +
geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
#geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
geom_label(aes(x=c(1:T), y=0, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
theme_classic(base_size=22) +
labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")
## Warning: Removed 31845 rows containing missing values or values outside the scale range
## (`geom_line()`).
ggplot()+
geom_jitter(data=df, aes(x=time, y=totalcount, group=simnumber), shape=".", width=0.3, height=0.3) +
geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
#geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
geom_label(aes(x=c(1:T), y=-0.5, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
theme_classic(base_size=22)+
labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")
## Warning: Removed 31845 rows containing missing values or values outside the scale range
## (`geom_point()`).
# ggplot()+
# geom_hex(data=df, aes(x=time, y=totalcount, group=simnumber), bins=10) +
# geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
# #geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
# geom_label(aes(x=c(1:T), y=-0.5, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
# scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
# theme_classic(base_size=22)+
# labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")
# ~~~~ Plot graph with population growth rates (Fig. 3.15) ~~~~
# op <- par(mar=c(4, 4, 2, 1), las=1, cex=1.1, "mfrow")
# layout(matrix(c(1, 1, 2, 3), 2, 2, byrow=TRUE), widths=c(1, 1), heights=c(1, 1), TRUE)
# plot(r[,1], type="l", lwd=0.5, ylab="Annual population growth rate", xlab="Time",
# ylim=range(r[which(!is.na(alive))]), col="lightgrey", axes=FALSE)
# axis(1); axis(2)
# for (s in 2:nsim){
# lines(r[!is.na(alive[,s]),s], lwd=0.5, col="lightgrey")
# }
# lines(apply(r, 1, mean, na.rm=TRUE), lwd=1.5)
# mtext("A", at=0, cex=1.5)
# a <- hist(mean.r, nclass=50, col="dodgerblue", main="",
# xlab="Population growth rate", axes=FALSE) #xlim=c(-2.5, 0.1),
# axis(1)
# axis(2, at=c(0, 5000, 10000, 15000, 20000, 25000, 30000), labels=c(0, 5, 10, 15, 20, 25, 30))
# mtext("B", at=a$mids[1], cex=1.5)
# a <- hist(mean.r[not.extinct], nclass=25, col="dodgerblue", main="",
# xlab="Population growth rate", axes=FALSE)
# axis(1)
# axis(2, at=c(0, 2000, 4000, 6000, 8000, 10000), labels=c(0, 2, 4, 6, 8, 10))
# mtext("C", at=a$mids[1], cex=1.5)
# par(op)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# https://github.com/mikemeredith/IPM_code/blob/main/IPM_03/IPM_03.3.5.R
# https://github.com/mikemeredith/IPM_code/blob/main/IPM_03/IPM_03.3.5.R
# 3.3 Classical analysis of a matrix population model
# ===================================================
# 3.3.5 Analysis of a matrix population model with different sources of
# stochasticity and parameter uncertainty
# ---------------------------------------------------------------------
# Define mean, measurement error and temporal variability of the demographic parameters
mean.sj <- 0.1163615 # Mean value of juv. survival
sd.sj.e <- 0.0554588 # Uncertainty of mean juv. survival as SD on natural scale
sd.sj.t <- plogis(0.03304038) # population variability of juv. survival as SD on logit scale # Andrew, this might have to be plogis(0.03304038), but not sure, check the text.
mean.sa <- 0.5222 # Mean value of ad. survival
sd.sa.e <- 0.05769693 # Uncertainty of mean ad. survival as SD on natural scale
sd.sa.t <- plogis(0.0386221) # Temporal variability of ad. survival as SD on logit scale
mean.f1 <- 1.605 # Mean value of productivity of 1y females
sd.f1.e <- 0.1326032 # Uncertainty of mean productivity as SD on natural scale
sd.f1.t <- 0.1014479 # Temporal variability of productivity as SD on natural scale
mean.fa <- 1.605 # Mean value of productivity of adult females
sd.fa.e <- 0.1326032 # Uncertainty of mean productivity as SD on natural scale
sd.fa.t <- 0.101 # Temporal variability of productivity as SD on natural scale
mean.imm <- 0.192587 # Mean value of adult immigration
sd.imm.e <- 0.06674809 # Uncertainty of mean adult immigration as SD on natural scale
sd.imm.t <- 0.02872281 # population variability of adult immigration as SD on natrual scale
# Define the number of years with predictions and the Monte Carlo setting
T <- 10 # Number of years (projection time frame)
nsim <- 10000 # Number of replicate populations simulated
# Define population matrix and initial stage-specific population sizes
N <- array(NA, dim=c(2, T+1, nsim))
N[,1,] <- c(3,1)
r <- matrix(NA, nrow=T, ncol=nsim)
alive <- matrix(NA, nrow=T, ncol=nsim)
mean.r <- numeric(nsim)
# Project population
for (s in 1:nsim){ # Loop over replicate populations
#if(s %% 100 == 0) {cat(paste("*** Simrep", s, "***\n")) } # Counter
# Generate a mean of the demographic rates (subject to measurement error)
msj <- rbeta2(1, mean.sj, sd.sj.e)
msa <- rbeta2(1, mean.sa, sd.sa.e)
mf1 <- rnorm(1, mean.f1, sd.f1.e)
mfa <- rnorm(1, mean.fa, sd.fa.e)
mia <- rnorm(1, mean.imm, sd.imm.e)
# Generate annual demographic rates (subject to temporal variability)
sj <- plogis(rnorm(T, qlogis(msj), sd.sj.t))
sa <- plogis(rnorm(T, qlogis(msa), sd.sa.t))
f1 <- pmax(0, rnorm(T, mf1, sd.f1.t)) # Avoids negative values
fa <- pmax(0, rnorm(T, mfa, sd.fa.t)) # Ditto
ia <- pmax(0, rnorm(T, mia, sd.imm.t)) # Ditto
# Project population (include demographic stochasticity)
for (t in 1:T){ # Loop over years
N[1,t+1,s] <- rpois(1, sj[t] * (f1[t] * N[1,t,s] + fa[t] * N[2,t,s]))
N[2,t+1,s] <- rbinom(1, sum(N[,t,s]), sa[t])+rpois(1,ia)
if (sum(N[,t+1,s]) == 0) break
#r[t,s] <- log(sum(N[,t+1,s])) - log(sum(N[,t,s]))
alive[t,s] <- t
} #t
#mean.r[s] <- mean(r[min(alive[,s], na.rm=TRUE):max(alive[,s], na.rm=TRUE),s])
}
#mean(mean.r)
#sd(mean.r)
not.extinct <- which(!is.na(alive[T,]))
mean(mean.r[not.extinct])
## [1] 0
sd(mean.r[not.extinct])
## [1] 0
# Extinction probability (after T years)
sum(is.na(alive[T,])) / nsim
## [1] 0.7736
rowSums(is.na(alive)) / nsim
## [1] 0.0300 0.1266 0.2507 0.3666 0.4616 0.5496 0.6284 0.6869 0.7358 0.7736
df <- data.frame(
totalcount=c(colSums(N[,,1:nsim])),
time=rep(c(0:T),nsim),
simnumber=gl(nsim, T+1)
)
# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# boxplot(c(colSums(N[,,1:nsim]))~factor(rep(c(0:T),nsim)))
bxpobj<-boxplot(c(colSums(N[,,1:nsim]))~factor(rep(c(0:T),nsim)), plot=F)
# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# points(bxpobj$stats[3,]~c(0:T), type="l", col="red", lwd=2)
#
# par(mai=c(1,1,0.1,1))
# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# points(bxpobj$stats[3,]~c(0:T), type="b", col="blue", lwd=2, pch=16)
# par(new=T)
# plot(c(rowSums(is.na(alive))/nsim)~c(1:T), type="b", col="red", bty="n", xlab="", ylab="", axes=F, ylim=c(0,1), xlim=c(0,T), lwd=2, pch=16)
# axis(side=4, at=pretty(range(c(0,1))))
# mtext("probability of extinction", side=4, line=3)
par(mai=c(1,1,0.1,1))
plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3, xlab="years of simulation", ylab="estimated population size")
points(bxpobj$stats[3,]~c(0:T), type="b", col="blue", lwd=2, pch=16)
par(new=T)
plot(c(rowSums(is.na(alive))/nsim)~c(1:T), type="b", col="red", bty="n", xlab="", ylab="", axes=F, ylim=c(0,1), xlim=c(0,T), lwd=2, pch=16)
axis(side=4, at=pretty(range(c(0,1))))
mtext("proportion of simulations reached N=0", side=4, line=3, col="red")
ggplot()+
geom_line(data=df, aes(x=time, y=totalcount, group=simnumber), linewidth=1, color="gray", alpha=0.8) +
geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
#geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
geom_label(aes(x=c(1:T), y=0, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
theme_classic(base_size=22) +
labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")
## Warning: Removed 38362 rows containing missing values or values outside the scale range
## (`geom_line()`).
ggplot()+
geom_jitter(data=df, aes(x=time, y=totalcount, group=simnumber), shape=".", width=0.3, height=0.3) +
geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
#geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
geom_label(aes(x=c(1:T), y=-0.5, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
theme_classic(base_size=22)+
labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")
## Warning: Removed 38362 rows containing missing values or values outside the scale range
## (`geom_point()`).
# ggplot()+
# geom_hex(data=df, aes(x=time, y=totalcount, group=simnumber), bins=10) +
# geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
# #geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
# geom_label(aes(x=c(1:T), y=-0.5, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
# scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
# theme_classic(base_size=22)+
# labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")
# ~~~~ Plot graph with population growth rates (Fig. 3.15) ~~~~
# op <- par(mar=c(4, 4, 2, 1), las=1, cex=1.1, "mfrow")
# layout(matrix(c(1, 1, 2, 3), 2, 2, byrow=TRUE), widths=c(1, 1), heights=c(1, 1), TRUE)
# plot(r[,1], type="l", lwd=0.5, ylab="Annual population growth rate", xlab="Time",
# ylim=range(r[which(!is.na(alive))]), col="lightgrey", axes=FALSE)
# axis(1); axis(2)
# for (s in 2:nsim){
# lines(r[!is.na(alive[,s]),s], lwd=0.5, col="lightgrey")
# }
# lines(apply(r, 1, mean, na.rm=TRUE), lwd=1.5)
# mtext("A", at=0, cex=1.5)
# a <- hist(mean.r, nclass=50, col="dodgerblue", main="",
# xlab="Population growth rate", axes=FALSE) #xlim=c(-2.5, 0.1),
# axis(1)
# axis(2, at=c(0, 5000, 10000, 15000, 20000, 25000, 30000), labels=c(0, 5, 10, 15, 20, 25, 30))
# mtext("B", at=a$mids[1], cex=1.5)
# a <- hist(mean.r[not.extinct], nclass=25, col="dodgerblue", main="",
# xlab="Population growth rate", axes=FALSE)
# axis(1)
# axis(2, at=c(0, 2000, 4000, 6000, 8000, 10000), labels=c(0, 2, 4, 6, 8, 10))
# mtext("C", at=a$mids[1], cex=1.5)
# par(op)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# https://github.com/mikemeredith/IPM_code/blob/main/IPM_03/IPM_03.3.5.R
# https://github.com/mikemeredith/IPM_code/blob/main/IPM_03/IPM_03.3.5.R
# 3.3 Classical analysis of a matrix population model
# ===================================================
# 3.3.5 Analysis of a matrix population model with different sources of
# stochasticity and parameter uncertainty
# ---------------------------------------------------------------------
# Define mean, measurement error and temporal variability of the demographic parameters
mean.sj <- 0.1163615 # Mean value of juv. survival
sd.sj.e <- 0.0554588 # Uncertainty of mean juv. survival as SD on natural scale
sd.sj.t <- plogis(0.03304038) # population variability of juv. survival as SD on logit scale # Andrew, this might have to be plogis(0.03304038), but not sure, check the text.
mean.sa <- 0.5222 # Mean value of ad. survival
sd.sa.e <- 0.05769693 # Uncertainty of mean ad. survival as SD on natural scale
sd.sa.t <- plogis(0.0386221) # Temporal variability of ad. survival as SD on logit scale
mean.f1 <- 0.343 # Mean value of productivity of 1y females
sd.f1.e <- 0.1326032 # Uncertainty of mean productivity as SD on natural scale
sd.f1.t <- 0.1014479 # Temporal variability of productivity as SD on natural scale
mean.fa <- 0.343 # Mean value of productivity of adult females
sd.fa.e <- 0.1326032 # Uncertainty of mean productivity as SD on natural scale
sd.fa.t <- 0.101 # Temporal variability of productivity as SD on natural scale
mean.imm <- 0.192587 # Mean value of adult immigration
sd.imm.e <- 0.06674809 # Uncertainty of mean adult immigration as SD on natural scale
sd.imm.t <- 0.02872281 # population variability of adult immigration as SD on natrual scale
# Define the number of years with predictions and the Monte Carlo setting
T <- 10 # Number of years (projection time frame)
nsim <- 10000 # Number of replicate populations simulated
# Define population matrix and initial stage-specific population sizes
N <- array(NA, dim=c(2, T+1, nsim))
N[,1,] <- c(3,1)
r <- matrix(NA, nrow=T, ncol=nsim)
alive <- matrix(NA, nrow=T, ncol=nsim)
mean.r <- numeric(nsim)
# Project population
for (s in 1:nsim){ # Loop over replicate populations
#if(s %% 100 == 0) {cat(paste("*** Simrep", s, "***\n")) } # Counter
# Generate a mean of the demographic rates (subject to measurement error)
msj <- rbeta2(1, mean.sj, sd.sj.e)
msa <- rbeta2(1, mean.sa, sd.sa.e)
mf1 <- rnorm(1, mean.f1, sd.f1.e)
mfa <- rnorm(1, mean.fa, sd.fa.e)
mia <- rnorm(1, mean.imm, sd.imm.e)
# Generate annual demographic rates (subject to temporal variability)
sj <- plogis(rnorm(T, qlogis(msj), sd.sj.t))
sa <- plogis(rnorm(T, qlogis(msa), sd.sa.t))
f1 <- pmax(0, rnorm(T, mf1, sd.f1.t)) # Avoids negative values
fa <- pmax(0, rnorm(T, mfa, sd.fa.t)) # Ditto
ia <- pmax(0, rnorm(T, mia, sd.imm.t)) # Ditto
# Project population (include demographic stochasticity)
for (t in 1:T){ # Loop over years
N[1,t+1,s] <- rpois(1, sj[t] * (f1[t] * N[1,t,s] + fa[t] * N[2,t,s]))
N[2,t+1,s] <- rbinom(1, sum(N[,t,s]), sa[t])+rpois(1,ia)
if (sum(N[,t+1,s]) == 0) break
#r[t,s] <- log(sum(N[,t+1,s])) - log(sum(N[,t,s]))
alive[t,s] <- t
} #t
#mean.r[s] <- mean(r[min(alive[,s], na.rm=TRUE):max(alive[,s], na.rm=TRUE),s])
}
#mean(mean.r)
#sd(mean.r)
not.extinct <- which(!is.na(alive[T,]))
mean(mean.r[not.extinct])
## [1] 0
sd(mean.r[not.extinct])
## [1] 0
# Extinction probability (after T years)
sum(is.na(alive[T,])) / nsim
## [1] 0.9372
rowSums(is.na(alive)) / nsim
## [1] 0.0554 0.2181 0.4013 0.5635 0.6808 0.7702 0.8315 0.8801 0.9137 0.9372
df <- data.frame(
totalcount=c(colSums(N[,,1:nsim])),
time=rep(c(0:T),nsim),
simnumber=gl(nsim, T+1)
)
# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# boxplot(c(colSums(N[,,1:nsim]))~factor(rep(c(0:T),nsim)))
bxpobj<-boxplot(c(colSums(N[,,1:nsim]))~factor(rep(c(0:T),nsim)), plot=F)
# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# points(bxpobj$stats[3,]~c(0:T), type="l", col="red", lwd=2)
#
# par(mai=c(1,1,0.1,1))
# plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3)
# points(bxpobj$stats[3,]~c(0:T), type="b", col="blue", lwd=2, pch=16)
# par(new=T)
# plot(c(rowSums(is.na(alive))/nsim)~c(1:T), type="b", col="red", bty="n", xlab="", ylab="", axes=F, ylim=c(0,1), xlim=c(0,T), lwd=2, pch=16)
# axis(side=4, at=pretty(range(c(0,1))))
# mtext("probability of extinction", side=4, line=3)
par(mai=c(1,1,0.1,1))
plot(jitter(c(colSums(N[,,1:nsim])),1)~jitter(rep(c(0:T),nsim),1), pch=3, xlab="years of simulation", ylab="estimated population size")
points(bxpobj$stats[3,]~c(0:T), type="b", col="blue", lwd=2, pch=16)
par(new=T)
plot(c(rowSums(is.na(alive))/nsim)~c(1:T), type="b", col="red", bty="n", xlab="", ylab="", axes=F, ylim=c(0,1), xlim=c(0,T), lwd=2, pch=16)
axis(side=4, at=pretty(range(c(0,1))))
mtext("proportion of simulations reached N=0", side=4, line=3, col="red")
ggplot()+
geom_line(data=df, aes(x=time, y=totalcount, group=simnumber), linewidth=1, color="gray", alpha=0.8) +
geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
#geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
geom_label(aes(x=c(1:T), y=0, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
theme_classic(base_size=22) +
labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")
## Warning: Removed 53146 rows containing missing values or values outside the scale range
## (`geom_line()`).
ggplot()+
geom_jitter(data=df, aes(x=time, y=totalcount, group=simnumber), shape=".", width=0.3, height=0.3) +
geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
#geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
geom_label(aes(x=c(1:T), y=-0.5, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
theme_classic(base_size=22)+
labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")
## Warning: Removed 53146 rows containing missing values or values outside the scale range
## (`geom_point()`).
# ggplot()+
# geom_hex(data=df, aes(x=time, y=totalcount, group=simnumber), bins=10) +
# geom_line(aes(x=c(0:T), y=bxpobj$stats[3,]), color="blue", linewidth=3)+
# #geom_line(aes(x=c(1:T), y=c(rowSums(is.na(alive))/nsim)), col="red", linewidth=2)+
# geom_label(aes(x=c(1:T), y=-0.5, label=round(c(rowSums(is.na(alive))/nsim), 2)), color="red")+
# scale_x_continuous(breaks=c(1:10))+scale_y_continuous(breaks=seq(0,50, by=2))+
# theme_classic(base_size=22)+
# labs(x="Years of simulated population size and\nproportion of simulations with N=0", y="Total simulated population count")
# ~~~~ Plot graph with population growth rates (Fig. 3.15) ~~~~
# op <- par(mar=c(4, 4, 2, 1), las=1, cex=1.1, "mfrow")
# layout(matrix(c(1, 1, 2, 3), 2, 2, byrow=TRUE), widths=c(1, 1), heights=c(1, 1), TRUE)
# plot(r[,1], type="l", lwd=0.5, ylab="Annual population growth rate", xlab="Time",
# ylim=range(r[which(!is.na(alive))]), col="lightgrey", axes=FALSE)
# axis(1); axis(2)
# for (s in 2:nsim){
# lines(r[!is.na(alive[,s]),s], lwd=0.5, col="lightgrey")
# }
# lines(apply(r, 1, mean, na.rm=TRUE), lwd=1.5)
# mtext("A", at=0, cex=1.5)
# a <- hist(mean.r, nclass=50, col="dodgerblue", main="",
# xlab="Population growth rate", axes=FALSE) #xlim=c(-2.5, 0.1),
# axis(1)
# axis(2, at=c(0, 5000, 10000, 15000, 20000, 25000, 30000), labels=c(0, 5, 10, 15, 20, 25, 30))
# mtext("B", at=a$mids[1], cex=1.5)
# a <- hist(mean.r[not.extinct], nclass=25, col="dodgerblue", main="",
# xlab="Population growth rate", axes=FALSE)
# axis(1)
# axis(2, at=c(0, 2000, 4000, 6000, 8000, 10000), labels=c(0, 2, 4, 6, 8, 10))
# mtext("C", at=a$mids[1], cex=1.5)
# par(op)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# https://github.com/mikemeredith/IPM_code/blob/main/IPM_03/IPM_03.3.5.R
# # for actual treatment fecundity values
# contm2 <- out1$mean$lambda[1] #2.575
# lowm2 <- out1$mean$lambda[7] #2.264
# medm2 <- out1$mean$lambda[14] #1.605
# highm2 <- out1$mean$lambda[20] #0.343
# alpha <- out1$mean$alpha #0.9422234
# beta <- out1$mean$beta #-1.073293
probofextinctdf <- data.frame(
treatment=c(rep("Control",10),rep("Low",10),rep("Medium", 10), rep("High",10)),
time=rep(1:10, 4),
probextinct=c(0.0241,0.0896, 0.1795, 0.2620, 0.3443, 0.4128, 0.4762, 0.5264, 0.5677, 0.6027,
0.0251, 0.1022, 0.2020, 0.2966, 0.3807, 0.4602, 0.5222, 0.5747, 0.6169, 0.6544,
0.0331, 0.1322, 0.2553, 0.3724, 0.4682, 0.5516, 0.6191, 0.6773, 0.7268, 0.7654,
0.0543, 0.2188, 0.4035, 0.5623, 0.6821, 0.7708, 0.8324, 0.8764, 0.9078, 0.9308),
fecundity=c(rep(2.575, 10),rep(2.264,10),rep(1.605,10),rep(0.343,10))
)
ggplot()+
geom_point(data=probofextinctdf, aes(x=time, y=probextinct, color=factor(treatment, levels=c("Control","Low","Medium","High"))))+
geom_line(data=probofextinctdf, aes(x=time, y=probextinct, color=factor(treatment, levels=c("Control","Low","Medium","High"))))+
labs(color="", x="Year of simulation", y="Proportion of simulations with N=0")+
scale_y_continuous(breaks=c(0, 0.25, 0.5, 0.75, 1.0), limits=c(0,1))+
scale_x_continuous(breaks=c(1:10))+
theme_classic()
plot(probextinct~fecundity, data=probofextinctdf)
abline(lm(probextinct~fecundity, data=probofextinctdf))
abline(lm(probextinct~fecundity, data=probofextinctdf|>filter(time==10)), col="red")
summary(lm(probextinct~fecundity+time, data=probofextinctdf))
##
## Call:
## lm(formula = probextinct ~ fecundity + time, data = probofextinctdf)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.20902 -0.02568 0.01506 0.03662 0.10939
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.226164 0.029695 7.616 4.35e-09 ***
## fecundity -0.123788 0.011829 -10.465 1.31e-12 ***
## time 0.079617 0.003527 22.572 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06407 on 37 degrees of freedom
## Multiple R-squared: 0.9436, Adjusted R-squared: 0.9406
## F-statistic: 309.5 on 2 and 37 DF, p-value: < 2.2e-16
summary(lm(probextinct~fecundity, data=probofextinctdf|>filter(time==10)))
##
## Call:
## lm(formula = probextinct ~ fecundity, data = filter(probofextinctdf,
## time == 10))
##
## Residuals:
## 1 2 3 4
## -0.0069874 -0.0008397 0.0136363 -0.0058092
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.986849 0.012874 76.66 0.00017 ***
## fecundity -0.146470 0.006773 -21.62 0.00213 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0116 on 2 degrees of freedom
## Multiple R-squared: 0.9957, Adjusted R-squared: 0.9936
## F-statistic: 467.6 on 1 and 2 DF, p-value: 0.002132
summary(lm(probextinct~fecundity, data=probofextinctdf|>filter(time==5)))
##
## Call:
## lm(formula = probextinct ~ fecundity, data = filter(probofextinctdf,
## time == 5))
##
## Residuals:
## 1 2 3 4
## 0.009575 -0.001511 -0.014634 0.006571
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.727902 0.014706 49.50 0.000408 ***
## fecundity -0.152690 0.007737 -19.73 0.002558 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01325 on 2 degrees of freedom
## Multiple R-squared: 0.9949, Adjusted R-squared: 0.9923
## F-statistic: 389.4 on 1 and 2 DF, p-value: 0.002558
# -0.20*(1/-0.1526)
# = 1.310616
# A reduction in fecundity of 1.31 is expected to lead to a 0.2 increase in risk of N=0 in the 5th of 10 years
# 2.5750723-1.310616=1.264456
# alpha = 0.9422234
# beta = -1.073293
# (log(1.264456)-0.9422234)/-1.073293 =
# 0.6592621